Thymic selection of T-cell receptors as an extreme value problem 
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T lymphocytes (T cells) orchestrate adaptive immune responses upon activation. T cell activation 
requires sufficiently strong binding of T cell receptors (TCRs) on their surface to short peptides (p) 
derived from foreign proteins, which are bound to major histocompatibility (MHC) gene products 
(displayed on antigen presenting cells). A diverse and self-tolerant T cell repertoire is selected 
in the thymus. We map thymic selection processes to an extreme value problem and provide an 
analytic expression for the amino acid compositions of selected TCRs (which enable its recognition 
functions) . 

PACS numbers: 87.10.-e, 02.50.-r,87.19.xw, 87.14.ep, 87.14.ef 



The adaptive immune system clears pathogens from 
infected hosts with the aid of T lymphocytes (T cells). 
Foreign (antigenic) and self-proteins are processed into 
short peptides (p) inside antigen-presenting cells (APC), 
bound to MHC proteins, and presented on the surface of 
APCs. Each T-cell receptor (TCR) has a conserved re- 
gion participating in the signaling functions, and a highly 
variable segment responsible for antigen recognition. Be- 
cause variable regions are generated by stochastic rear- 
rangement of the relevant genes, most T cells express 
a distinct TCR. The diversity of the T cell repertoire 
enables the immune system to recognize many different 
antigenic short pMHC complexes. Peptides presented 
on MHC class I are typically 8-11 amino acids long [l|, 
which is enough to cover all possible self-pcptides (the 
human proteome consists of P « 10 7 amino- acids @, Q) 
as well as many antigenic peptides. TCR recognition 
of pMHC is both specific and degenerate. It is specific, 
because most mutations to the recognized peptide amino 
acids abrogate recognition 0, [f| • R is degenerate because 
a given TCR can recognize several antigenic peptides Q . 

The gene rearrangement process ensuring the diversity 
of TCR is random. R may thus result in T cells poten- 
tially harmful to the host, because they bind strongly 
to self peptide-MHC complexes; or useless T cells which 
bind too weakly to MHC to recognize antigenic pep- 
tides. Such aberrant TCRs are eliminated in the thy- 
mus 0, [H, 0, [n| , where immature T cells (thymocytes) 
are exposed to a large set (10 3 — 10 4 ) of sclf-pMHC. Thy- 
mocytes expressing a TCR that binds with high affinity 
to any sclf-pMHC molecule are deleted in the thymus 
(a process called negative selection). However, a thy- 
mocyte's TCR must also bind sufficiently strongly to at 
least one self pMHC complex to receive survival signals 
and emerge from the thymus (a process called positive 
selection). 

Signaling events, gene transcription programs, and 
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13, 14, 15, lal have been studied exten- 



sively. Despite many important advances, how interac- 
tions with self-pMHC complexes in the thymus shape the 
peptide-binding properties of selected TCR amino acid 
sequences, such that mature T cells exhibit their special 
properties, is poorly understood. To address this issue, 
in Ref. 17[ we numerically studied a simple model where 



TCRs and pMHC were represented by strings of amino 
acids (Fig. Q]) . These strings indicate the amino-acids on 
the interface between TCRs and pMHC complexes, and 
it is assumed that each site on a TCR interacts only with 
a corresponding site on pMHC. The binding interface of 
TCR is actually composed of a region that is in contact 
with the MHC molecule, and a segment that is in con- 
tact with the peptide. It is the latter part that is highly 
variable, while the former is more conserved. We shall 
therefore explicitly consider only the former amino-acids, 
but not the latter. Similarly, there are many possible 
peptides that can bind to MHC, and their sequences are 
considered explicitly, whereas those of the MHC are not 
(there are only a few types of MHC in each individual 
human [l[). We could in principal add a few sites to the 
TCR and pMHC strings to account for any variability in 
the segments not considered. 

Simplified representations of amino-acids (e.g., as a 
string of numbers or bits) were employed earlier [l5|, EH 
[lli | in the context of TCR-pMHC interactions, mainly 
to report that negative selection reduces TCR cross- 
reactivity. In Ref. [17[ , we numerically studied the model 
in Fig. [T] (and described below) to qualitatively describe 
the role of positive and negative selection on the amino- 
acid composition of selected TCRs. By randomly gen- 
erating TCR and pMHC sequences, and implementing 
thymic selection in silico, we showed that selected TCRs 
are enriched in weakly interacting amino acids, and ex- 
plained how this leads to specific, yet cross-reactive, TCR 
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FIG. 1: Schematic representation of the interface between 
TCR and pMHC complexes. The segment of TCR that is 
in contact with peptides is highly variable and modeled by a 
string of N amino-acids. The peptide is also modeled by a 
sequence of length N, and the binding energy is computed as 
a sum of pairwise interactions. We don't explicitly consider 
TCR sites in contact with MHC, as they are more or less 
conserved, and only assign them a net interaction energy E c . 

recognition of antigen, a long-standing puzzle. In this 
paper we show that the model can be solved exactly in 
the limit of long TCR/peptide sequences. The result- 
ing analytic expression for the amino-acid composition 
of selected TCRs is surprisingly accurate even for short 
peptides and provides a theoretical basis for previous nu- 
merical results. Furthermore, we are able to obtain a 
phase diagram that indicates the ranges of parameters 
where negative or positive selection are dominant, lead- 
ing to quite different bias in selection/function. 

To assess the effects of thymic selection, as well as anti- 
gen recognition, we evaluate the free energy of interaction 
between TCR-pMHC pairs (for brevity, free energy will 
be referred to as energy). The interaction energy is com- 
posed of two parts: a TCR interaction with MHC, and a 
TCR interaction with the peptide. The former is given a 
value E c (which may be varied to describe different TCRs 
and MHCs). The latter is obtained by aligning the TCR 
and pMHC amino-acids that are treated explicitly, and 
adding the pairwise interactions between corresponding 
pairs. For a given TCR-pMHC pair, this gives 

N 

E int (t,s)=E c + J2J(ti,Si), (1) 

?:=i 

where J(ti,Si) is the contribution from the ith amino 
acids of the TCR (U) and the peptide (sj), and N is the 
length of the variable TCR/peptide region. The matrix J 
encodes the interaction energies between specific pairs of 
amino-acids. For numerical implementations we use the 
Miyazawa-Jernigan (MJ) matrix that was developed 
in the context of protein folding. 

Immature T cells interact with a set S of M sclf-pMHC 
complexes, where typically M is of the order of 10 3 — 10 4 . 
To mimic thymic selection, sequences that bind to any 
self-pMHC too strongly (E int < E n ) are deleted (nega- 
tive selection). However, a thymocyte's TCR must also 
bind sufficiently strongly (i?i n t < E p ) to at least one 
sclf-pMHC to receive survival signals and emerge from 



the thymus (positive selection). A thymocyte expressing 
TCR with string t will thus be selected if the strongest in- 
teraction with sclf-pMHC is between thresholds for neg- 
ative and positive selection, i.e. 

E n <mm{E int (t,s)} < E p . (2) 

Recent experiments [ll| show that the difference between 
thresholds for positive and negative selection is relatively 
small (a few ksT). 

Equati on j[2l) casts thymic selection as an extreme value 
problem [20(, enabling us to calculate the probability 
P se \(t) that a TCR sequence t will be selected in the 
thymus. Let us indicate by p{x\t) the probability den- 
sity function (PDF) of the interaction energy between the 
TCR f and a random peptide. The PDF Il(x\t) of the 
strongest (minimum) of the M independent random in- 
teraction energies is then obtained by multiplying p with 
the probability of all remaining (M — 1) energy values be- 
ing larger- (l — P[E < x\t)} 1 , where P(E < x\t) is 
the cumulative probability- and noting the multiplicity 
M for which energy is the lowest. The probability that 
TCR t is selected is then obtained by integrating H{x\t ) 
over the allowed range, as 

Pm\(t) = [ "ll(x\t)dx, With 

Je„ 

n(x\t) = M p(x\t)(l-P(E <x\t)) M ^ . (3) 

For M>1, this extreme value distribution (EVD) con- 
verges to one of three possible forms, (20j depending on 
the tail of the PDF for each entry. Equation (QJ indicates 
that in our case as each energy is the sum of N contribu- 
tions, p{x\t) should be a Gaussian for large N, in which 
case the relevant EVD is the Gumbel distribution. [2(| 

To obtain an explicit form for n(a;|t), we model the 
set S of self-pcptidcs as M strings in which each amino- 
acid is chosen independently. The probability f a for 
selecting amino-acid a at each site is taken to be the 
frequency of this amino-acid in the self-proteomc. For 
a specific TCR sequence t, the average interaction en- 
ergy with self peptides follows from Eq. (JT|) as E av (t ) = 
E c + Z)i=i £(*»)> with £(ti) = [J(t t ,a)} a , where we 
have denoted the average over self amino-acid frequen- 
cies by [G(a)] a = Y^ a =i faG(a). Similarly, the variance 
of the interaction energy is V(t) = ^ =1 V(ti), where 
V(U) = [J(t i: a) 2 ] a - [J(U,a)] 2 a . For large N, we can 
approximate p{x\t ) with a Gaussian PDF with the above 
mean and variance. From standard results for the Gum- 
bel distribution [2(j, we conclude that in the limit of 
M>1, the peak of the distribution H(x\i ) is located at 

E (t) = E &v {t) - ^2V{t)\nM , (4) 

and its width is E (t) = y / 7r 2 V r (f)/(121nM). (Since the 
PDF p(x\t ) originates from a bounded set of energies, 
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it is strictly not Gaussian in the tails. Hence, once the 
extreme values begin to probe the tail of the distribu- 
tion, the above results will no longer be valid. Indeed, in 
the limit when M ~ 0(20^), the EVD will approach a 
delta-function centered at the TW-independent value cor- 
responding to the optimal binding energy.) 

In the limit of long TCR/pcptidcs (TV 1), we can 
exactly calculate the statistics of the amino-acid com- 
position of selected TCRs. To obtain a proper thermo- 
dynamic limit, we need to set {E c , E p , E n } cx TV, and 
hxM cx TV. The latter ensures that the peak of the dis- 
tribution, Eo(t), is proportional to TV, and also results 
in a width T> (t) which is independent of TV. (The re- 
lation In A/ = aTV can be justified with the expectation 
that M should grow proportionately to the protcomc size 
P, while TV oc In P to enable encoding the proteome.) In 
this large TV limit, the EVD is sufficiently narrow that the 
value of the optimal energy can be precisely equated with 
the peak Eo(t), and Eq. ([2|) for the selection condition 
can be replaced with 



corresponds to a set of non-interacting variables, with 
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21nM^V0i) < E p . (5) 



Thus, for each sequence t, we have to evaluate the 'Hamil- 
tonian' E n (t), and the sequence is accepted if this energy 
falls in the interval (E n ,E p ). This is somewhat similar 
to the micro-canonical ensemble in Statistical Physics, 
with the restriction of the energy to an interval rather 
than a narrow range only a minor elaboration (see be- 
low). From the equivalence of canonical and micro- 
canonical ensembles for large TV, we know that the proba- 
bility for a sequence is governed by the Boltzmann weight 
p(t) cx (Hill /ti) exp[-/3E (t)]. Here {f a } indicate the 
natural frequencies of the different amino-acids prior to 
selection, while the effect of thymic selection is captured 
in the parameter (3 which is determined by solving for 
the average energy. 

The appearance of y / 2\nMJ2i V(U) in the Hamilto- 
nian initially appears as a complication that makes exact 
computation of the average energy from exp[— (3Eo(t )] 
impossible. However, this apparent 'coupling' is easily 
dealt with by standard methods such as Legendre trans- 
forms or Hamiltonian minimization [211 ]. This can be 
justified easily as follows: We need to solve a 'Hamilto- 
nian' 7i(U, V) which depends on two extensive quanti- 
ties U = and V = Ej!LiV(t<). The corre- 
sponding partition function can be decomposed as Z = 
v ri{U,V)e~P n ^ U ' V \ but can be approximated with 
its largest term. Note that the same density of states 
fi(t/, V) = e S{u,v)/k B a pp earS; irrespective of the specific 
form of H(U, V). In particular, the choice Ho — E c + U — 
7 F-lnM/(2 7 ) = £ c +Eti[(£&)-7V&)]-lnM/(2 7 ) 
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p(t) cx Y[ fu exp 
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(3(£(U) ~ 7 V(*0) 



(6) 



for which thermodynamic quantities (such as entropy) 
are easily computed. By judicious choice of 7 we can 
then ensure that the same average energy appears for 
Tto(t) and our E Q (t). Using Legendre transforms, which 
is equivalent to minimizing Ho(t) with respect to 7 , 
one finds that the required Eq (t ) is obtained by setting 

7 (/3) = y / lnM/(2TV (V)^), where (• • • )^ refers to the 

average with the non-interacting weight e~^ u ~ lV ' . 

Finally, the value of (3 has to be determined by con- 
straining the average energy determined above to the 
range in Eq. ((5]), while maximizing entropy. Given the 
bounded set of energies, the inverse temperature (3 can be 
either negative or positive. The 20 N possible values for 
Eo(t ) span a range from £" m i n to E max , and a correspond- 
ing number of states which is a bell-shape between these 
extremes with a maximum at some E m id- If ^mid > E p , 
we must set j3 such that (Eo(t)\ = E p . In this case, 
f3 > 0, positive selection is dominant and stronger amino- 
acids are selected. If TS m id < E n , we must set (3 such that 
(Eo(t )) = E n , p < 0, negative selection is dominant and 
weaker amino-acids are selected. For E n < £"mid < E p , 
we must set (3 = and there is no modification due to 
selection. 

Figure [2] depicts the variation of (3 as a function of 
ln(M)/TV and threshold for negative selection E n with 
(E p - E n )/N = 0.5k B T. Consider TCRs that do not 
bind too strongly or weakly to MHC, as such TCRs are 
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FIG. 2: (Color) Color representation of the dependence of the 
inverse temperature (3 on the number of self-peptides In M/N 
and the threshold for negative selection energy E n /N with 
(E p - E„)/N = 0.5k B T in the limit of large TV. The region 
between the black lines corresponds to j3 = 0, to the right 
(left) of which negative (positive) selection is dominant, and 
weak (strong) amino-acids are selected. Note that as (E p — 
E n )/N goes to zero, the intermediate region disappears. The 
dotted lines indicate the relevant parameter values for thymic 
selection in mouse (see text) that result in (3 = — 0.37(fcsT) _1 . 
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(captured by the parameter 7) which was not elucidated 
from the limited numerical data. Furthermore, the phase 
diagram in Fig. [2] indicates how upon raising the num- 
ber of self-pcptides there is a transition from preference 
for strong amino-acids (/? > 0, positive selection dom- 
inant) to weak amino-acids (/? < 0, negative selection 
dominant), which may be feasibly tested in future exper- 
iments, along the lines in Ref. [4|. 

This work was supported by National Institutes of 
Health (NIH) Grant 1-PO1-AI071195-01 and a NIH Di- 
rector's Pioneer award (to A.K.C.). 



FIG. 3: Amino-acid composition of selected TCR sequences, 
ordered in increasing frequency along the abscissa. The data 
points in black are obtained numerically with the parameters 
relevant to mouse (see caption of Fig. [2] and text). The error 
bars reflect the sample size used to generate the histograms 
and differences for different realizations of M self-peptides. 
The dashed line is the result of the EVD analysis in the large 
N limit from Eq. (J7J, and the agreement is quite good. In both 
cases we have used the Miyazawa-Jernigan matrix J H p. an d 
amino acid frequencies f a from the mouse proteome 0, [l7[ . 



unlikely to be selected (e.g., E n - E c = -21k B T). For 
the set of parameters that are relevant for thymic selec- 
i.e. N = 5 

-0.37(fc s T) 



tion in mouse 
M = 10 3 , we find /3 



E p — E, 



= 2.hk B T and 
which means that 



negative selection is dominant and weaker amino-acids 
are selected. Also 7 = 0.83(ksT)~ 1 indicating a pref- 
erence for amino-acids with smaller variations in bind- 
ing energy. With these parameters we can calculate the 
amino-acid frequencies of selected TCRs as 



•I a 



(sel) 



/ a exp[-^(g( g )- 7 V( g ))] 
E£i A «p [-/?(£(&) - 7 V(6))] 



(7) 



It is important to ask if the above expression, exact in 
the limit of iV — > 00, has any relevance to the actual 
TCR/peptides with N ~ 5 — 10. We thus numerically 
simulated the case of N = 5 by generating a random 
set of 10 6 TCR sequences, and selected them against 
M = 10 3 self-pcptides. The selected TCRs were used to 
construct the amino-acid frequencies depicted in Fig. [3] 
The dashed line in this figure comes from Eq. |(7J) with 
the same J and {f a }- The agreement between the two is 
remarkable given the small value of N = 5, and may be 
indicative of small corrections to the N — > 00 result. 

Equation (|7|) thus provides an analytical expression 
that captures the characteristics of TCR amino-acids se- 
lected against many peptides in the thymus. In accord 
with previous numerical results [l7j], and some available 
data from normal mouse, and human, it predicts (since 
f3 < 0) that TCR sequences are enriched in weakly inter- 
acting amino-acids (small £). This result was used pre- 
viously to explain their specificity. However, Eq. (J7]) 
further indicates the role of promiscuity of amino-acids 
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